Glow discharge in low pressure plasma PVD: mathematical model and 

numerical simulations 



A. Speranza, A. Monti 
Industrial Innovation Through Technological Transfer - PT^ Onlus 
alessandro . speranza@i2t3 . unif i . it 
L. Barletti, I. Borsi, L. Meacci 
Dip. di Matematica "U. Dini" Universita degli Studi di Firenze 

S. Fanfani 
Galileo Vacuum Systems srl, Prato 

December 14, 2010 



Abstract 

In this paper we analyze the problem of glow discharge 
in low pressure plasma in industrial plant, for chambers 
of different shapes and various working parameters, like 
pressure and electric potential. The model described 
is based upon a static approximation of the AC con- 
figuration with two electrodes and a drift diffusion ap- 
proximation for the current density of positive ions and 
electrons. A detailed discussion of the boundary con- 
ditions imposed is given, as well as the full description 
of the mathematical model. 

Numerical simulations were performed for a simple 
ID model and two different 2D models, corresponding 
to two different settings of the industrial plant. The 
simpler case consists of a radially symmetric chamber, 
with one central electrode (cathode), based upon a DC 
generator. In this case, the steel chamber acts as the 
anode. The second model concerns a two dimensional 
horizontal cut of the most common plant configuration, 
with two electrodes connected to an AC generator. The 
case is treated in a "quasi-static" approximation. The 
three models show some common behaviours, partic- 
ularly including the main expected features, such as 
dark spaces, glow regions and a wide "plasma region" . 
Furthermore, the three shown models show some simi- 
larities with previously published results concerning ID 
and simplified 2D models, as well as with some prelim- 
inary results of the full 3D case. 

keyv^rords: PVD, glow discharge, surface coating, 
plasma, drift-diffusion 

1 The industrial problem 

Plasma PVD (Physical Vapour Deposition) is a widely 
used process, applied for, particularly, in the coat- 



ing industry. The process basically aims at forming a 
thin coating layer on the surface of different substrates. 
PVD is normally used as an alternative physical pro- 
cess to classical wet chemical process, such as galvanic 
baths. 

The case we particularly refer to concerns the met- 
alization of plastic headlights for automotive (Fig. [1]) , 
with a thin aluminium layer, and its coating with thin 
polymeric film, such as HMDSO (Hexa Methyldisilox- 
ane). The HMDSO is used to protect the metalization 
from external agents and ageing. 

As it can be seen from Fig. [l]-b, the industrial plant 
consists of three main parts. The two semi-chambers 
visible at the bottom of the figure, contain the head- 
lights to be processed. The two moving parts close 
alternatively on the fixed part of the plant (rectangular 
in Fig. [T]-b) as shown in the photograph (Fig. [T]-a). In 
this configuration, the two moving semi-chambers con- 
tain circular plates disposed all around the axis of sym- 
metry, moving in planetary motion. The headlights are 
deposed on the plates during the process. The fixed 
part of the plant, visible as rectangle at the top of 
Fig. [IJ-b, contain the electrodes (see below) and the 
vacuum pump with all the pipelining. The industrial 
plant is normally made in two different configurations. 
The configuration represented in figure, has two cylin- 
drical vertical electrodes and works on alternate cur- 
rent (AC), i.e., an alternate electric field is applied to 
the electrodes in order to produce the glow discharge. 
In the second configuration, one cylindrical electrode 
(normally the cathode) is located in the center of the 
chamber and a direct electric potential (DC) is applied 
to it. In this second configuration, the metallic chamber 
itself works as the anode. 

The industrial process proceeds through three main 
steps 
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Figure 1: a) Photograph of the industrial plant by Galileo Vacuum Systems s.p.a., objecto of the work. The headlights 
on the left hand side, have been already exposed to metalization, while the semi-chamber on the right hand side is 
closed and the process is active, b) Schematic representation of the PVD chamber (not the one produced by Galileo 
Vacuum Systems) as seen from the top. The chamber has two semi-chambers which can contain the headlights to 
treat, deposed on circular plates. The rectangular part above the two semi-chambers contains the electrodes and the 
vacuum pump. 



• Vacuum: one of the two mobile semi-chambers is 
locked against the fixed part, and vacuum (~ 10"'^ 
mbar) is made 

• Metalization: aluminium filaments located all 
along the height of the mobile semi-chamber are 
heated via Joule effect. In such a low pressure at- 
mosphere, aluminium vaporizes and a thin layer 
of metallic aluminium deposits onto the headlights 
surface. The more uniform and denser is the cloud 
of vaporized aluminium, the more uniform and 
thicker will be the metallic layer on the surface. 

• Polymer coating: a strong electric potential differ- 
ence 2-5 kV) is applied to the electrodes. In 
these conditions, the residual atmosphere trans- 
forms into plasma and a glow discharge appears. 
At the same time, a monomer flows through noz- 
zles within the chamber, increasing the pressure 
up to l-5xl0~^ mbar. The plasma atmosphere in- 
duces the polymerization on the substrate of the 
headlights, and a thin polymer film covers the alu- 
minium layer. 

In this work, we will neglect completely the first two 
steps of the process, and will concentrate on the third 
one. In particular, we will concentrate on the part in 
which, in low pressure atmosphere, the electric poten- 
tial is applied to the electrodes (or to the cathode, in 
the configuration with one central electrode) and the 
glow discharge is induced [6l HH [HI EH [22l [28] , while 
the monomer flows and polymerization occurs. The aim 



of our work is to model the glow discharge process, in 
order to gain information concerning, e.g., the distri- 
bution of positive and negative charges, or the electric 
fleld within the chamber, during the glow discharge. It 
is clear, in fact, that HMDSO polymerizes under the 
influence of the plasma atmosphere. However, it is not 
clear in what particular conditions the polymerization 
is actually favored. In other words, our final aim is to 
describe the plasma atmosphere all across the chamber, 
in order to correlate its features, with the main features 
of the resulting polymer film, such as thickness, homo- 
geneity etc. With this information, one could try to 
optimize the shape of the chamber and the value of 
the main process parameters, in order to have the best 
possible result, in the widest part of the chamber. 

2 The physical and mathemati- 
cal model 

As mentioned in the previous section, this work focuses 
on the final step of the industrial process, i.e., on the 
polymerization of the HMDSO in low pressure plasma. 
In fact, our model concerns more the plasma atmo- 
sphere in which the HMDSO polymerizes, rather than 
the process of polymerization itself. The polymeriza- 
tion would be rather difficult to model, as the process 
is strongly influenced by different factors, not least the 
substrate on which the polymer film is formed or its 
shape, besides the conditions of the plasma itself (see. 



for instance, dZHIl]). For all these reasons, given that 
a realistic model of polymerization would be rather dif- 
ficult to achieve, and even more complex to solve nu- 
merically (due to the strongly asymmetric geometry of 
the domain), we focus on the analysis of the conditions 
in which the monomer polymerizes while flowing across 
the chamber. Our final aim is to describe the conditions 
of the plasma, in terms of electric field, density distri- 
bution of positive and negative charges etc., in order to 
correlate these conditions with the main features of the 
obtained polymer film. 

2.1 Hypotheses of the model 

Our mathematical model is based upon a set of simple 
assumptions, aimed at simplifying the set of equations 
and, in parallel, the geometry of the domain. As far as 
the chamber is concerned, we therefore assume 

• The steel chamber is perfectly cylindrical. 

• The electrodes (one or two, depending on the con- 
figuration of the plant) are perfectly cylindrical 
and not connected anyhow to the chamber. In 
other words, the electrodes "float" in the chamber, 
as they are generally shorter than the chamber it- 
self. 

In order to simplify the set up of the mathematical 
model, we then make the following assumptions, con- 
cerning the plasma atmosphere: 

• Thermal (and thermodynamic) equilibrium is as- 
sumed throughout the chamber. In this way, we 
assume a perfectly static atmosphere. 

• Plasma is made only of two species of charged par- 
ticles: positive single ionized atoms (ions) and free 
electrons. Clearly this assumption is rather strong, 
especially in the final part of the process, i.e., when 
the monomer is flowing within the chamber and the 
polymer film is forming. When this happens, the 
atmosphere is made of a great number of species, 
such as monomers, ions, radicals, pieces of polymer 
chains and so on, however, this is the only possible 
starting point for any kind of modelling approach. 

• We always consider a stationary electric field, even 
in the case of an alternate potential difference ap- 
plied to the electrodes. In this way we can neglect 
all the terms involving the variation of the electric 
field in the Maxwell equations (see below). 

Moreover, we will make a further assumption, concern- 
ing, as we will sec in later sections, the boundary condi- 
tions to impose to our set of equations: we assume that 
all the steel parts composing the chamber and the elec- 
trodes, are "perfectly absorbing", i.e., all the charged 



particles hitting the surface of any steel part, go back 
into the chamber as neutral particles. In other words, 
charged particles are never reflected into the chamber, 
as they are. 

Before proceeding to the description of the mathe- 
matical model, let us go back for a moment to the as- 
sumption of a static electric field mentioned earlier, in 
spite of the choice of the chamber configuration with 
one (DC generator) or two (AC generator) electrodes. 
In fact, this assumption is reasonable, as, even in the 
case of two electrodes, the frequency of the applied field 
is of the order of 100 kHz. Given that the radius of the 
chamber is normally of the order of 1 m, and that the 
wavelength of the applied field is of the order of 3 km, 
it turns out that, within the extension of the chamber, 
the electric field can be considered uniform. In other 
words, the chamber is sufficiently small to allow us the 
use of a "quasi-static approximation", i.e., to neglect 
all the possible delay effects due to the variation of the 
electric field within it. 

As we will shortly see in the following sections, in 
spite of the drastic simplifications of the physical prob- 
lem that we made in order to solve numerically the 
mathematical model, our results turn out to correlate 
fairly well with observations and direct measures. In 
fact, although the real atmosphere consists of many 
types of ions and radicals, resulting from the polymer- 
ization process, observations show, as in simpler atmo- 
sphere, the formation of a "dark space" surrounding 
the electrodes, confined by a bright region: "the glow" . 
As we will see in the following sections, in our numeri- 
cal results the two regions are clearly evident. The two 
regions differ in particular in the values of the number 
density of charged particles and electric field intensity. 
In other words, in spite of the high degree of approx- 
imation of our model, our numerical results turn out 
to give some reasonable insights of the actual problem, 
marking an important step forward in the interpreta- 
tion of observations and direct measures, aimed at the 
optimization of the industrial setting. 

2.2 Set up of the equations 

Once the physical model has been set as described in 
the previous section, the set of equations to solve turns 
out to be rather simple. As mentioned, we assume a 
steady atmosphere and thermal equilibrium across the 
chamber. Furthermore, we assume a steady electric 
field applied to the electrodes. In these hypotheses, 
the electric potential must obey the non homogeneous 
Poisson equation, i.e., 

Al/ = -^(n, -Tie) (1) 

where A is the Laplace operator, V is the electric po- 
tential, e is the electron electric charge (~ 1.6 10~^^ C), 



e is the absolute dielectric constant of the atmosphere 
and Ui, Tie are the number densities of ions and elec- 
trons respectively. For the equation above, one must 
set boundary conditions which, in this case, turn out 
to be straightforward Dirichlet boundary conditions, as 
it will be shown in the following section. 

As far as ions and electrons current densities Ji ^ 
are concerned, in the hypothesis of a continuous model 
(see below), they obviously must obey the stationary 
continuity equation, i.e., 

V-Ji = ^ieneS{V) (2) 
V-Je = ^ieneS{V) (3) 

where ji^ is the electrons mobility coefficient and S{V) 
is the ionization frequency, which is obtained from the 
Townsend formula [11 |35l [39] 

5(F)=Acxp|-b(^^^ |p|VF| (4) 

The term on the right hand side of Eqs. (j2l3p . 
fien-eS{V), basically expresses the number density of 
ion/electron pairs formed through a collision between 
a neutral particle and a free electron, per unit time. 
As obvious, this term depends on the pressure P and 
the electric field module |Vy|. Note that, compared 
to the standard notation {e.g., [HUHDi we pulled the 
electron mobility out of the definition of the Townsend 
coefficient, in order to make the typical scaling of the 
equations more evident. The two coefficients A and B 
in the equation above, are obtained by a fit of exper- 
imental data, and depend on the composition of the 
atmosphere (see, e. 3., [Ml [Ml 1301 133 [32]) • Typical val- 
ues of the two coefficients range between 4 and 60 cm~^ 
Torr~^, for A, and between 14 and 36 (V/cm Torr)°-^, 
for B. In the following, we made a fairly conservative 
choice, in order to control more easily the numerical 
calculations and used A between 2 and 3 cm~^ Torr~^ 
and B between 12 and 13 (V/cm Torr)" "', comparable 
with the values for He reported by Ward [32] and with 
the values used by Boeuf and Pitchoford [5] . As shown 
in Fig. [21 the first Townsend coefficient is strongly de- 
pendent on the choice of the set of coefficients A and B. 
However, as the electric field grows, and S{V) / P\W\ 
settles on a constant value, the frequency of ionization 
S{V)^e becomes linear with |VT^|. 

The current density Jq of charged particles q (i or e) 
in Eqs. (|2l3p is modeled through a simple drift-diffusion 
equation, i.e., 

Jq ;= -sgn{q)^ignqW - DgVug (5) 

where Dq is the diffusivity coefficient and is linked to 
the charge carrier mobility via the so called Einstein's 



formula 



where Kg is the Boltzmann constant, and Tq is the tem- 
perature of the charged particle. In the following, we 
made the assumption that ions are in thermal equilib- 
rium with the neutral gas, therefore Ti ^ T = 300 K, 
while the electron temperature was set to Te = 10000 
K throughout our domain, corresponding to an average 
energy Ee ^ 0.4 cV. The electrons and ions mobility 
are taken, coherently with the choice we made for the 
first Townsend coefficient, as fie — 1333 m^ V^^s^^ 
(/ieP=106 cm^V-^s-iTorr) and =10.66 m^ V-^s-^ 
(/^iP=8xlO'^ cm^V~^s~-^Torr) respectively, at a pres- 
sure of 10 Pa, as reported by Ward [55] . 

Our set of equations therefore consists of Eqs. ([T])- 
([2]) and ([3]) , with the use of ([5]) , in the three unknowns 
V, rii and Ue. 

2.3 Comments on the model 

As shown, our model is based upon the heuristic as- 
sumption, to be verified a posteriori, that a fluid model 
is suitable throughout the whole domain. Clearly, this 
assumption will have to be checked against the param- 
eters used to define the conditions of our plasma. It is, 
in fact, well known, that in a very low pressure plasma 
a collision-free sheath could be formed in contact with 
an absorbing or a neutral wall [T] [34] , effectively con- 
fining the neutral plasma. As the hypothesis of high 
collisionality is not satisfied, a continuous model may 
not be applicable to this case. Within the sheath re- 
gion, electrons arc almost absent and rii ^ rig. Fur- 
thermore, the sheath is confined by a transition region 
(the pre-sheath) where the plasma is almost exactly 
neutral (i.e., rii — n^). The sheath is stable when the 
so called Bohm criterion is satisfied, i. e. , when the ions 
drift velocity is larger than some critical speed, corre- 
sponding, as shown later [H [1] , to a Mach surface for 
the ions flux. As mentioned, when the conditions for 
the formation of the sheath are satisfied, a purely con- 
tinuous model cannot be used, as the assumption of 
high collisionality is not satisfied within the sheath and 
the pre-sheat. In this case, one possible approach is to 
develop a fiuid model that connects, in the limit, with 
the sheath model (see, 6.17., [121 [38]). However, the 
conditions for the formation of a stable sheath are not 
completely clear yet. It was shown that the Bohm cri- 
terion results, in fact, in a sufficient but not necessary 
condition for the sheath formation [37] [35] . It was also 
shown by Valentini [35] and Riemann [37], that such a 
criterion has to be strongly modified when collisions are 
present. In other words, a sheath could be formed even 
in coUisional plasma, although its characteristics are 
quite different in this case. Godyak and Sternberg [T3] 
showed that the model upon which the formation of 
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Figure 2: First Townsend coefficient rescaled by pressure and electron mobility as defined in Eq. ([4]), as a function 
of |VF|. The two sets of coefficients A ^ 2 — 3 cm^^ Torr^^ and i? = 12 — 13 (V/cm Torr)" "' correspond roughly to 
the data published in [39] for He. 



such a sheath is described, is intrinsically inconsistent, 
suggesting that a modification on the basic assump- 
tions should be made. Valentini |37j . on the other hand, 
showed that in a discharge plasma, i.e., a slightly ion- 
ized, low pressure plasma, similar to ours, the sheath 
can still be formed, under appropriate conditions. How- 
ever, the whole sheath concept, in particular in front of 
the cathode, should be somehow redefined. In fact, 
in these conditions, a fluid model can usually be used 
throughout the whole domain, as collisions are frequent 
and important, even though, the electrode fall still re- 
gion turns out to show very peculiar characteristics, 
compared to the neutral plasma region. In order to take 
into account of the peculiar behaviours of that region, 
authors chose to use different methods. Hagelaar et al., 
for instance |18j , used boundary conditions to take into 
account of the different behaviour of the fast electrons 
generated within the cathode fall region. Surendra et 
al. [33] . coupled a fluid model (similar to ours) with an- 
other model for the fast electrons. On the other hand, 
other authors llll[ni[Il[I3[ini[lil[lSll3Illin] used a 
fluid model throughout the whole domain, investigating 
its range of applicability, and choosing boundary con- 
ditions in order to have a tractable and yet physically 
meaningful model. Boeuf and Pitchford [5] argued that 
a fluid drift-diffusion model can normally be used with 
values of the pressure of the order of 100 mTorr 13 
Pa). In fact, the same approach was also used for RF 
glow discharge models [3J UJ (TSl [HI [131 [H] , as long as 
the number of collisions remains high enough, i.e., as 
long as the applied pressure is sufficiently high. The 



threshold value for the pressure normally grows from 
100 to 500 mTorr (~ 67 Pa), for applications of the 
fiuid model to an RF discharge model. However, the 
argument, in both cases, the AC and DC, is based upon 
the fact that, for pressures above some values, the ions 
mean free path is normally the smallest length scale, 
compared to, for instance, the electrons Debye length, 
the typical length scale of the model, and so on. In 
these conditions, collisions can be considered frequent 
enough to allow for a continuous model to be physically 
meaningful. 

A similar argument applies exactly also to our case. 
As mentioned above, we chose a pressure value of 10 Pa 
75 mTorr), which is very close to the lowest value 
mentioned, for instance by Boeuf and Pitchford [5] . In 
these conditions, the ions mean free path turns out to 
be of the order of 1 mm. This value should be compared 
with the typical length scale of the chamber (L = 1 
m) and, for instance, the electrons Debye length Aq. 
In our case, the electrons Debye length turns out to 

be Aj) = L / \J Te/fiej where Te = T^/T and hf. is the 
non dimensional electrons number density (see next sec- 
tion) . Since we set Te ~ 30 and from our numerical re- 
sults, we have at most hf. ~ 10^, we get that Ap > 10 
mm. In any case, near the electrodes, where fi^ decays 
down to 10'^, we have Aj) ^ 0.1 m. This particular 
result indicates [5^ that the electric field may not be 
strong enough to confine the ions in a very narrow re- 
gion in contact with the electrode. Instead, collisions 
allow the ions to diffuse further inside the bulk region. 
This general consideration is supported by our numer- 



ical results. As we will see in the following sections, 
the cathode sheath region, shows a dominance of ions, 
as it should be, but a still large presence of electrons 
(rie > 10^ m^'^). Furthermore, our sheath is confined 
by a region with a high net charge density, and not 
by a neutral region such as the one assumed to con- 
fine the sheath region in the collision-less plasma. In 
fact, as already predicted by Valentini [37], in our case, 
the net charge density reaches a maximum before de- 
caying rapidly to zero in the neutral plasma region, 
coherently with the assumption of a highly collisional 
plasma. Sheridan and Goree [32] showed that the non- 
coUisional and the collisional regimes connect smoothly 
through an intermeidate regime with the relavant index 
being, once again, the collision parameter a = / . 
The authors argued that the collisional regime extends 
in the region a > 0.1. For values of a of the order 10, 
like in our case, the regime is already highly collisional 
and the fluid model should be fully acceptable. 

In spite of the low pressure, we therefore conclude 
that our plasma shows a still strongly collisional regime 
that allows us, although close to the limiting case, to 
use a continuous model throughout the domain. How- 
ever, even when a smooth connection between the the 
numerical results of the fluid model and an asymptotic 
analysis of the boundary layer was attempted |38| . the 
two solutions showed an appreciable difference only at 
a very short length-scale which are far beyond our aims. 
Clearly, the argument of high collisionality applies only 
since we chose to focus on the final step of the industrial 
process (that we referred to as "Polymer coating" in 
Sec.[T]). The pressure there grows rapidly between 10^'^ 
mbar to 0.1-0.5 mbar as the monomer flows through the 
chamber and polymerization occurs. The same argu- 
ment, on the other hand, would not apply if we wanted 
to describe the conditions of the plasma atmosphere at 
the very start of the process, i.e., when the pressure is 
still 100 times lower, and the atmosphere consists only 
of the residual air present in the chamber. 

3 The problem of boundary con- 
ditions 

To the set of equations described in the previous sec- 
tion, one has to impose boundary conditions on each 
one of the steel parts which confine the problem do- 
main, i.e., the two electrodes and the chamber sur- 
face. As mentioned in the previous section, to the Pois- 
son equation for the electric potential Eq. ([T]), one can 
simply impose Dirichlet boundary conditions. For in- 
stance, defining Vc, the potential of the steel chamber, 
as the reference potential, and redefining all the poten- 
tials as potential differences V = V — Vc and swapping 



V and V, we can set 

V = Vq on the anode 

V = — Vb on the cathode 

V ~ on the steel chamber (6) 

In this way, the potential difference between the two 
electrodes is set to be 2Vb. The potential difference 
with the chamber is, respectively, Vq and — Vb for the 
anode and the cathode. 

There is no such a simple way to express boundary 
conditions for the number densities rii and rig, unless 
we make some further assumptions on the actual inter- 
action between charged particles and the electrodes or 
the chamber. Let us therefore analyze the condition to 
impose on the different parts of the plant, for the two 
species of charged particles. For the ions following, e.g., 
the arguments by Hagelaar et al. [18] , one can make the 
following assumptions: 

• The cathode absorbs ions, in the sense that any ion 
hitting the cathode surface, captures a free elec- 
tron emerging on the surface and is reflected into 
the chamber as a neutral particle. We assume, 
as mentioned earlier, that the cathode is a "per- 
fect absorber of ions" in the sense that none of the 
ions hitting the surface is simply reflected into the 
chamber. In these conditions, the density flux Jj 
reduces to just the drift term, while the diffusion 
vanishes, i.e., we set a simple Neumann boundary 
condition such as 

Vn^-v^O (7) 

where v is the normal versor pointing inward the 
cathode surface (i.e., outwards from our domain) 

• The anode does not emit (or absorb) ions, thus the 
net inward flux across its surface must vanish. In 
other words, on the anode we set a mixed boundary 
condition of the type 

- Ji • 1/ := -{-n^n,VV - D^Vui) • i/ = (8) 

where u is the normal versor pointing inward the 
anode surface (i.e., outwards from our domain) 

• Following this argument, the steel chamber will 
have to behave in a way which is similar to the 
anode or the cathode, depending on the direction 
of the electric field (and thus, of the ions flux). 
When the electric fleld —W "enters" the surface, 
i.e., has the same sign as the outwards normal ver- 
sor, the steel chamber behaves like a cathode. On 
the other hand, when the electric field —W is out- 
ward the surface, i.e., entering our domain, and 
therefore anti-parallel to the outward normal ver- 



sor, the steel chamber behaves hke the anode. In 
compact form, we can write this condition as 

Ji 1/ := {-fiiUiVV - DiVrii) ■ v = 

= —ai^iUiW ■ u (9) 

where = 1 (cathodc-Uke) if —W ■ v > Q and 
fli = (anode hkc) otherwise 

A symmetric argument appHes for the electrons. In 
this case, however, the roles of anode and cathode are 
inverted, thus 

• The anode absorbs electrons, in the sense that 
any electron, hitting the anode, is captured by it. 
Thus, the diffusion term of the electrons flux van- 
ishes on its surface, i.e., 

VUe • IV = (10) 

• On the cathode surface, the electrons inward flux 
is due to thermal and secondary emission, thus 

-Je-IV := -(/-te^eVy - DeVrie) • f = 

= 7»Ji ■ ^"e-^th (11) 

where 7^ is the secondary emission coefficient (0.1 
in our case) and w^j^ = y^Sn-gT/irm is the thermal 
velocity and m is the electron mass. 

• As in the case of the ions, the steel chamber be- 
haves differently, depending on the sign of the elec- 
tric field, i.e., 

Je • IV ;= (/^erie VV^ - D^VUe) ■ V = 

= ae/ie"eVy • V (12) 

where Oe = (cathode- like) if — VV^ • iv > and 
ae ~ I (anode like) otherwise 

3.1 Final setting of the model 

In order to make our result as general as possible, let 
us define new non dimensional variables u ~ V/Vq for 
the potential, hi = Ui/no and fie = Ue/nQ, where Vq 
and no are scaling factors to be chosen appropriately. 
Let us also rescale all the coordinates by the typical 
length scale of the chamber (L = 1 m). Inserting the 
new variables into Eqs. (IT]|3]) and swapping hi and hg 
with Hi and rie, respectively, we can set Vq = k-qT /e ('^ 
lO-^V) and no = en^T/e^L^ (- lO^m'^). With this 
notation, we get our final model in non dimensional 
form as 

— Au = Ui — rig (13) 
V • (-n,Vu - Vui) = —neSiu) (14) 



V • ^ngVit — TeVne j = neS(u) (15) 
with boundary conditions, in compact form 

— [—sgn{q)nqVu — V?ig] • v = 
= sgn(q) a n^Vu ■ v + 6ge (^iJi ■ f + ^"g^'th^ (1^) 

where we have defined a non dimensional electron tem- 
perature Tp = T^/T, and we simply set S{u) = S{Vqu). 
In the above equations, q indicates the charge carrier, 
i.e., i or e, and a = 1 if — sgn(g)Vu ■ u > and a = 
otherwise. The last term starting with the Kronecker 
delta 6qe applies only to the electrons on "cathode-like" 
boundary conditions, and is neglected elsewhere. 

4 Numerical results 

Let us now show some numerical results obtained for 
some simple configurations of the chamber in ID and 
2D. All the simulations were performed with Comsol® 
Multiphysics [3], by using a modified convection- 
diffusion model, for the two current densities, coupled 
with an electrostatic model (Poisson equation) for the 
electric potential. 

4.1 One dimensional model 

In order to be able to interpret correctly the results 
obtained in 2D, particularly in the configuration with 
two electrodes, let us start with showing some results 
obtained in ID. In this case, our domain is simply a 
straight line extending between the cathode, on the 
left, and the anode, on the right. The steel chamber 
in this case has no particular role and, consequently, 
the boundary conditions greatly simplify as there is no 
influence of direction of the electric field on the surface. 
In spite of this drastic reduction of the model, as we 
will see, the results already contain an important part 
of the information given by the 2D case, in more realis- 
tic configurations. Let us then start with a plot of the 
electric potential and electric field norm. In Fig. |3] we 
represent the non dimensional electric potential u (left) 
and its spatial derivative u^, i.e., the electric field (with 
a positive sign), in logarithmic scale. The electric po- 
tential clearly shows two main regimes of the solution. 
In the first regime corresponding to a region extend- 
ing between the cathode and the center of the domain, 
the electric potential grows steeply. Next to this, in a 
region extending from the center to the other end of 
the domain, at the anode, the electric potential has al- 
ready reached its boundary value. However, when we 
plot the electric field (right), in logarithmic scale, we 
actually observe a very rich behaviour, showing all the 
main expected features. First of all, starting from the 
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Figure 3: Non dimensional electric potential (left) and norm of the electric field (right). The dimensional potential 
difference between anode and cathode is set at 5000 V. The norm of the electric field shows clearly the main features 
expected, such as the "cathode dark space" , the negative and positive glow and a plasma region, with low electric 
field. 



cathode (left), the "cathode dark space" extends for 
roughly one fourth of the domain. The dark space is a 
region relatively poor of ions (compared to the plasma 
region), where the charged particles are being acceler- 
ated by a strong electric field. In the chamber, this 
small region appears clearly as a dark space extending 
all along the cathode. Starting from the right end of 
the dark space, up to almost a half of the domain, we 
can observe the "negative glow" . In this region, a large 
number of highly energetic ions (accelerated in the dark 
space by the electric field), release their energy in form 
of light. In the experimental setting, one would see a 
vivid purple zone all around the cathode. Next to the 
negative glow is a "plasma region" . In this region, al- 
most nothing happens, in the sense that, in spite of 
the large number of charged particles (see Fig. |4]) , the 
electric field here is so small that the energy density 
turns out to be very low. Furthermore, the number 
of ions and electrons perfectly balance, thus making 
the plasma region practically neutral (while the cath- 
ode dark space and the negative glow are positively 
charged. Fig. |3]-right). The other end of the domain, 
confined by the anode, is occupied by the positive glow 
(and the anode dark space). In this region the electric 
field increases again, although not as much as in the 
cathode dark space. However, as shown in Fig. HI in 
this region, the electrons exceed the ions, making the 
plasma negatively charged. 

The representation of the net charge of the plasma, 
in Fig. |4]-right confirms all our assertions. The cathode 
dark space, is positively charged, and extends from the 
cathode up to roughly 0.25. In the negative glow region. 



the net charge reaches its maximum and then decreases 
steeply towards the plasma region, which is neutral as 
a consequence of the perfect overlapping of the number 
densities of ions and electrons. Finally, the positive 
glow shows a negatively charged atmosphere, as the 
number of electrons exceeds the number of ions in the 
nearby of the anode. Overall, Fig. U shows that the 
plasma is positively charged as the total number of ions 
exceeds the total number of electrons, except in the 
nearby of the anode. 

4.2 Two dimensional model — 1 elec- 
trode 

Let us now show some partial numerical results in the 
two dimensional model. In order to correlate our results 
to the previous ID case, let us assume a simple con- 
figuration with one central electrode (cathode), while 
the steel chamber acts like the anode and is kept at 
a neutral potential. The cathode is set at a potential 
(dimensional V) of -3000 V with respect to the steel 
chamber. 

In this simple case (Fig. [5]), we can use the cylindrical 
symmetry of the domain, and reduce the real 3D case 
to a simpler 2D model. In fact, one can reduce the inte- 
gration domain to just a half (we chose the lower half) 
of the vertical section of the chamber, as the upper half 
is perfectly symmetric to the lower part. The follow- 
ing results were obtained with a FEM triangular mesh 
(Fig. [5]) consisting of 37456 non uniform elements. The 
system of equations was solved through an affine invari- 
ant form of the damped Newton method [9l [10] , imple- 




Figure 4: Left: non dimensional number density of ions and electrons (dashed). In the insert, same plot in logarithmic 
scale. As already shown in the plot of the electric field, the main expected features such as the "cathode dark space" , 
the "negative glow", and the "positive glow" are evident. In the dark space, the number of ions is greatly larger 
than the number of electrons, but not yet maximum. A maximum is reached in the glow region. The plasma region 
is neutral, as the number of ions balances exactly the number of electrons (overlapping curves). In the positive glow 
region, the number of electrons exceeds the number of ions and the plasma in negatively charged. Right: net charge 
density p = ni — ne- The existence of the four different regions is evident as the net charge density changes drastically 
from one region to the other. 



mcnted through a direct linear solver UMFPACK [S] 
included in Comsol® Multiphysics. 

As in the previous case, a plot of the non dimensional 
electric potential (Fig. E]), shows the two main regimes 
with the potential increasing steeply from the cathode 
inside the region, and a wide plasma region, where the 
electric potential has already reached a value which is 
very close to the one imposed at the anode (u = in 
this case). The similarity between the behaviour of the 
ID case and the 2D case with one central electrode is 
evident when we plot the electric potential evaluated 
along a horizontal line drawn over the symmetry axis 
that closes the domain at the top (Fig. |6]-right). As in 
the previous ID model, the two regimes (actually hid- 
ing the four regions already shown in the ID model) 
with the potential increasing steeply proceeding from 
the cathode inwards the integration domain, and the 
"plasma region" with an almost constant potential, are 
clearly evident. A similar parallelism between the ID 
model and the cylindrically symmetric 2D one is evi- 
dent when we plot the number densities of ions Ui and 
electrons rig (Fig. [71). Again, as in the previous case, 
the two distributions overlap perfectly in the "plasma 
region" , while near the cathode the number of ions is 
largely greater than the number of electrons. Vicev- 
ersa, near the anode (the steel chamber in this case), 
the electrons exceed the ions, thus making the plasma 
locally negative, but positive overall. Fig. [7] shows also 
a clear resemblance with the number density of positive 
ions and electrons reported by Boeuf in 0. The main 



differences between our case and theirs, besides, pos- 
sibly, the value of the scaling factors, is due probably 
due to the different choice of boundary conditions for 
the potential of the steel chamber (linear in their case, 
uniform in ours). 

4.3 Two dimensional model — two elec- 
trodes 

Let us now approach the case with the widest indus- 
trial interest, as it resembles, although yet in a two 
dimensional representation, the real plant, in the con- 
figuration with two electrodes and an alternate electric 
field. As the cylindrical symmetry is lost, due to the 
introduction of the second electrode, none of the pos- 
sible two dimensional reductions can represent the real 
three dimensional model. However, a horizontal slice 
of the three dimensional domain (Fig. |8]), as we will 
see, already gives some useful insights of what hap- 
pens in the central (in height) area of the electrodes of 
the full three dimensional case, where we can assume 
some sort of limited translational invariance. The FEM 
mesh consisted of 6560 non uniform standard triangu- 
lar elements. As in the previous cases in one or two 
dimensions, the 2D model with two electrodes shows 
the typical expected behaviour. Once again, the plot of 
the electric potential (contour plot in Fig. IHl-left) shows 
a dominance of the region extending from the anode 
inwards, where the potential is roughly constant. A 
plot of the potential evaluated along the orthogonal 
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Figure 5: Representation of the two dimensional mesh used in the cylindrically symmetric case with one electrode. 
The domain is delimited on the left by the cathode, while the steel chamber (bottom and right) constitutes the 
anode. At the top we impose a symmetry condition as on the central axis. 



cross section between the two electrodes (sketched in 
Fig. [3]-lcft as a horizontal segment between the elec- 
trodes) shows that the two-dimensional case resembles 
accurately the previous ID and 2D models. As above, 
the electric potential remains constant in the nearby of 
the anode and decreases steeply only as we get closer to 
the cathode. The plot of the norm of the electric field 
(non dimensional) , in Fig. |9]-right /bottom, shows again 
the features already observed in the ID model, i.e., the 
dark space, the two glows and the plasma region, al- 
though the extension of the regions is now completely 
different from the previous cases. 

In spite of the orthogonality between the two projec- 
tions (horizontal in this case vs vertical in the cylindri- 
cally symmetric one), the great resemblance of this 2D 
case with the previous cylindrically symmetric model, 
is clearly evident also when we plot the number density 
of ions rii and (Fig. [TO)) . However, in this case, the 
number density shows some extreme features that were 
not so evident in the previous cases. Fig. [TOl-top shows 
that both the number density of ions (left) and elec- 
trons (right) are strongly peaked in the plasma region. 
In fact, as before, the two peaks are exactly overlap- 
ping, as it is evident from the cross section plot of the 
two densities (bottom- left) and of the net charge ni — Ue 
(bottom- right), where we see how, in the plasma re- 
gion, the net charge is, once again, almost exactly null. 
As before, the plasma is then positively charged in the 
cathode dark space and the cathode glow, and nega- 
tively charged in the anode glow. 

Finally, let us go back for a moment to the bound- 
ary conditions imposed in to the steel chamber in 
Eqs. (|9ll2p . There, we said, we have to impose a differ- 
ent condition, depending on the sign of the electric field 
on the surface. In other words, we assumed that, some- 



how, the steel chamber contributes to the total current 
crossing the chamber, behaving, in some parts, anode- 
like, and in others, cathode-like. An insight of this ef- 
fect, which is clearly rather weak, but yet present, can 
be obtained by showing a plot of the stream lines of the 
current density of ions J^. Fig. 1111 shows, as expected, 
that the ions flux (and similarly the electrons fiux, not 
shown) is concentrated within the region between the 
electrodes. However, some (non negligible) density is 
fiowing also in the area around the two electrodes, be- 
tween the electrodes and the steel chamber. There, the 
electric field streamlines close between the electrodes 
and the steel chamber. In other words, the steel cham- 
ber, through a partial coupling with the electrodes, con- 
tributes to the flow of electric current. In particular, 
the steel chamber behaves like an anode in the area 
where the streamlines connect it to the cathode, while 
it behaves like a cathode, where the streamlines con- 
nect it to the anode (most of the surface in this case). 
As mentioned, the contribution of this partial coupling 
to the total amount of current is fairly low. However, 
this contribution exists and confirms that the bound- 
ary conditions must include it somehow. Furthermore, 
this suggests that also a coupling of different nature, for 
instance, due to residual capacity effects, is possible be- 
tween the steel chamber and the electrodes. This type 
of coupling was completely neglected in our case, since 
we assumed a static field, however, this result shows 
that, in case of alternate current, it could contribute 
even further to the passage of electric current through 
the plant. 
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Figure 6: Left: Representation of the non dimensional electric potential u in the two-dimensional case with one 
central electrode. The potential increases steeply from the cathode into the domain, reaching quickly a value already 
close to the one imposed at the anode u — 0. The arrows pointing inwards the cathode show the strength of the 
electric field in that area, compared to the plasma region, where they are invisible as the electric field is low. Right: 
Cross section obtained by drawing a horizontal line at the top of the domain, i.e., on the horizontal symmetry axis. 
The plot shows the perfect correspondence of this case, to the ID case. 



5 Conclusions 

We analyzed the behaviour of plasma atmosphere in 
a PVD industrial plant by means of a mathematical 
model involving the electric potential and the number 
density of ions and electrons. The drastic simplifica- 
tions we made in setting the mathematical model, al- 
lowed us to solve it numerically in a fairly straightfor- 
ward way. The main difficulty of the model, as already 
mentioned in the literature, concerns the setting of the 
boundary conditions. In fact, previous results show 
that the usual Neumann or Dirichlet boundary condi- 
tions are too simple to give physically meaningful re- 
sults. For this reason, we made an intermediate choice 
of mixed boundary conditions, that take into account 
of, at least, secondary and thermal emission of elec- 
trons. However, as in our case the plant configuration 
consists of two electrodes and the steel chamber, the 
boundary conditions turn out to depend on the solution 
inside the domain. The three equations of the model 
are therefore coupled in a strongly non linear way. In 
spite of this, for fairly simple geometries, such as the 
shown two dimensional cut, and the corresponding full 
three dimensional case (not included), the numerical 
solution of the model is still possible. In fact, the two 
dimensional model with two electrodes show a great re- 
semblance with the one dimensional case and the cylin- 
drically symmetric two dimensional case. The three 
analyzed cases show a similar behaviour of the electric 
potential, with evidence of two different regimes. A 



plasma region is clearly evident with almost constant 
potential, bounded by the anode. A highly energetic 
glow region with the potential increasing steeply, ex- 
tends between the plasma and the cathode. 

Similar considerations can be made, concerning the 
number densities of ions and electrons. In all the shown 
cases (and similarly it happens in the three dimensional 
case, although that is rather more difficult to visualize) 
the number density of ions and electrons overlap almost 
exactly in the plasma region. In the cathode dark space 
and positive glow, the ions greatly dominate, while in 
the negative glow (close to the anode), the number of 
electrons exceeds the number of ions. This behaviour 
was again expected, since it was both observed exper- 
imentally and simulated by a mathematical model for 
the one dimensional case (such as the glass tube [5]) 
and some simple configurations of the 2D model [6]. 
As far as we know, this is the first case where a full 
two-dimensional model, with no rotation invariance, is 
solved numerically and fully reported. Furthermore, as 
far as we know, it is the first time that a clear compari- 
son between the ID model and more complex 2D mod- 
els, either with one or two electrodes, is made. Clearly, 
this is not conclusive, as the full three dimensional was 
not analyzed in detail. However, it is the first time 
that a published numerical model actually takes into 
the account of the real industrial configuration of an 
alternate current (AC) plant. 

In terms of the application of our results in an in- 
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Figure 7: Left: Contour plot (top) and 3D surface (bottom) of the non dimensional number density of ions (top). 
Right: Contour plot (top) and 3D surface (bottom) of the non dimensional number density of electrons. As in the ID 
model, the two distributions overlap perfectly in the central area (plasma region) of the domain. Near the cathode 
(dark space and glow) the number of ions largely exceeds the number of electrons, thus making the plasma overall, 
positively charged. 



dustrial contest, one should try to correlate these the- 
oretical considerations with some experimental results 
of the deposition process, in order to show the influ- 
ence of the different conditions of the plasma on the 
polymeric film formed. Clearly, only partial consider- 
ations can be made by a discussion of our numerical 
results. One expects, for instance, the plasma region to 
be fairly quiet, in terms of electric energy, as in spite 
of the great number of charged particles (both ions and 
electrons reach their absolute maximum in the plasma 
region), the electric field is almost null. Conversely, in 
the negative glow region, a large number of ions (com- 
pared to the electrons) are strongly accelerated by a 
strong electric field. This is the most energetic region 
of the whole domain, and one expects a rather dramatic 
behaviour of the plasma. For example, in the glow re- 
gion (or the cathode dark space), one expects a strong 



"sputtering effect" due to the ions hitting the surface 
of substrate to be treated and, possibly, removing the 
aluminium film. That is clearly a phenomenon that one 
would want to avoid as much as possible, when choos- 
ing where to pose the pieces to be treated within the 
chamber. However, at this stage, i.e., without a sys- 
tematic experimental campaign it is hard to give any 
further interpretation of the model, except for very gen- 
eral considerations, concerning the real industrial pro- 
cess of polymerization and coatings. In spite of this, 
we believe our model can be a good starting point for a 
systematic investigation of the real process. In order to 
interpret correctly our results and guide the industrial 
process, one should, possibly, not only compare them 
with experiments and direct observations, but also with 
the results of a simplified polymerization model. One 
could then correlate the results from the two models 
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Figure 8: Representation of the mesh used to solve the model in the two dimensional case, obtained by cutting 
horizontally the original 3D domain. The two cylindrical electrodes are represented by two circles, the anode on the 
left, and the cathode on the right. 



with observations, and draw some conclusions about 
the real process of coating. For instance, a deeper in- 
vestigation of the polymerization process could suggest 
that a highly energetic region is favorable in terms of 
speed and effectiveness of the polymerization. 

Although this project of investigation is just at the 
beginning, we believe that our model marked a good 
step forward in the direction of a more systematic and 
scientific approach to the optimization of a commonly 
used, but not yet fully understood, industrial process. 
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Figure 10: Non dimensional density of ions (top/left) and electrons (top/right) in three dimensional representation. 
The two distributions are strongly peaked in the plasma region. However, the two peaks overlap almost exactly, 
as clear from the cross section plot (bottom/left) and from the cross section plot of the net charge density rii — Ue 
(bottom/right). The net charge density (bottom/right) shows also the typical behaviour encountered in the previous 
cases. The plasma region is roughly neutral, while the cathode dark space and the negative glow are positively 
charged. The positive glow is negatively charged, as the number of electrons exceeds the number of ions. 
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Figure 11: Ions density flux Ji and streamlines of the electric field. The current flows, as expected, mainly between the 
two electrodes, with a high density of ions flux in the vicinity of the cathode. However, some electric field streamlines 
close up on the steel chamber. The steel chamber therefore contributes to the total current flow within the chamber, 
due to a partial coupling between the electrodes and the chamber itself. In the area where the streamlines connect 
the chamber to the anode, the chamber behaves like a cathode, while in the area where the streamlines connect to 
the cathode, the chamber behaves like the anode. 
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